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ABSTRACT 

We use high-resolution cosmological zoom-in simulations from the Feedback in Realistic En¬ 
vironment (FIRE) project to study the galaxy mass-metallicity relations (MZR) from z = 0-6. 
These simulations include explicit models of the multi-phase ISM, star formation, and stel¬ 
lar feedback. The simulations cover halo masses Mhaio = 10^-10^^ Mq and stellar masses 
M* = lO^-lO" Mq at z = 0 and have been shown to produce many observed galaxy prop¬ 
erties from z = 0-6. Eor the first time, our simulations agree reasonably well with the ob¬ 
served mass-metallicity relations at z = 0-3 for a broad range of galaxy masses. We pre¬ 
dict the evolution of the MZR from z = 0-6, as log(Zgas/Z 0 ) = 12-f log(0/H) — 9.0 = 
0.35 [log(M*/MQ) - 10] -f 0.93exp(-0.43z) - 1.05 and log(Z*/ZQ) = [Ee/H]-|-0.2 = 
0.40 [log(M*/M 0 ) — 10] -|-0.67exp(—0.50z) — 1.04, for gas-phase and stellar metallicity, re¬ 
spectively. Our simulations suggest that the evolution of MZR is associated with the evolution 
of stellar/gas mass fractions at different redshifts, indicating the existence of a universal metal¬ 
licity relation between stellar mass, gas mass, and metallicities. In our simulations, galaxies 
above = 10 ^ Mq are able to retain a large fraction of their metals inside the halo, because 
metal-rich winds fail to escape completely and are recycled into the galaxy. This resolves a 
long-standing discrepancy between “sub-grid” wind models (and semi-analytic models) and 
observations, where common sub-grid models cannot simultaneously reproduce the MZR and 
the stellar mass functions. 
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1 INTRODUCTION 

The galaxy mass-metallicity relation (MZR) is one of the most 
fundamental properties observed in galaxies. In the local universe, 
there is a tight correlation between galaxy stellar mass and gas- 
phase oxygen abundance for star-forming galaxies (e.g. |Tremonti| 
|et al. 12004] ), with an intrinsic scatter of only 0.1 dex in log(0/H). 
This relation has been extended to local dwarf galaxies and found 
to be a uniform, tight correlation over five orders of m agnitude in 
stellar mass, from M, = 10^-10*’ Mq < Lee et al. 20061. Many dif- 


[ 2008 |pa 


et al.|[2009^, despite the fact that the results are 


ferent groups have confirmed the MZR to exist not only in the local 
universe but also at high redshifts up to z ~ 2.3 (e.g. Savaglio et al. 


2005 

|Erb et al.|20061|Zahid et al.|20111|2012 

Andrews & Martini 

2013 

IHenry et al.|2013a|bHYabe et al.|2014 

Steidel et al.|20141 


Sanders et al.|20l5| l. Zahid et al. 1 20131 compiled a number of the 

observed MZR from z = 0-2.3 and found that the MZR evolves 
with redshift, with higher metallicity at low redshift for a given 
stellar mass. The MZR is also found at z > 3 (e.g. |Maiolino et al.| 


obtained from very small samples. 


Gas-phase metallicities represent the “current” state of chem¬ 
ical enrichment in the galaxies, while stellar metallicities reflect 
the “time-averaged” galactic metallicity across the whole star for¬ 
mation history. Similarly, an MZR is also found in stellar metal¬ 
licities. [Gallazzi et al.| P005^ derived the stellar metallicities for 
~44,000 galaxies from SDSS and found a tight correlation between 
stellar mass a nd stellar metallici ty for galaxies of stellar masses 
lO^-lO'^ Mq. Kirby et al. (2013 i measured the metallicities of in- 
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dividual stars in a sample of dwarf galaxies within the Local Group 
and found the SDSS stellar MZR can be continually extended down 
to 10^ Mq. Despite the fact that stellar metallicity is challenging to 
measure at high redshifts, the stellar MZR provides very important 
and complimentary insights on the chemical evolution of galaxies, 
especially for massive quiescent galaxies and satellite galaxies in 
the local group where the gas-phase metallicities are hard to mea¬ 
sure due to their low gas content. 

Simple analytic models of galactic chemical evolution, such 
as the “closed-box”, “leaky-box”, and “accreting-box” models (e.g. 
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Schmijl|]^3[ rfalbot & Amett|197r||Searle & Sargent|1972[[Ed-| 


munds|19^ i, are often quoted to illustrate the qualitative behavior 

of the MZR. More complicated models have also been developed 
to work in cosmological contexts and to connect gas inflows, out- 


flows, and star formation to galactic chemical evolution (e.g. Dal- 

canton|2007 

Einlator & Dave|20081 |Dave et al.|2012||Lilly et al. 

20131|Lu eta] 

.|2015a^. These models indicate that the existence of 


MZR is the consequence of an interplay between star formation ef¬ 
ficiency, metal loss from gas outflows, and gas recycling and accre¬ 
tion. For example, the stellar mass-halo mass relation (e.g. |Moster| 
|et al.|20l^ [Behroozi et al.|2013| > indicates that the star formation 
efficiency (fraction of baryons turned into stars) is lower in low- 
mass galaxies than in more massive galaxies, suggesting that low- 
mass galaxies should be less metal-enriched. Meanwhile, galactic 
winds are ubiquitous (see e.g. [Veilleux et al.|[2005| for a recent 
review), carrying metals away from galaxies. Low mass galaxies 
have shallow potential wells so they tend to lose a significant frac¬ 
tion of their gas and metals, while massive galaxies have potential 
wells deep enough to prevent material from escaping the galaxy 
(e.g. |Dekel & Silk|p'98^ . On the other hand, gas inflows bring 
the metal-poor gas in the galactic halo and/or in the intergalactic 
medium (IGM) inwards, diluting the metal content in the interstel¬ 
lar medium (ISM) and supplying new material for star formation 
(e.g. |Keres et al.|2Q05[ [Faucher-Giguere et al.|201l) . During this 
process, a considerable fraction of the gas and metals that have been 
formerly ejected via outflows eventually come back to the galaxy 
(e.g. |Bertone et al.|2007[|0ppenheimer et al.|2010| l. Galaxy merg¬ 
ers and AGN activity could also be important, in the sense that they 
can trigger violent starburst, drive intensive gas outflows, and ulti¬ 
mately quench the star formation in the galaxy (e.g.|Springel et al.| 
|2005|[Ho^ns et al.|2013'^ . 

Analytical models usually rely on simplified assumptions such 
as perfect mixing and adopt simple analytic prescriptions describ¬ 
ing star formation, gas accretion, and outflows. In reality, these 
physical processes are tightly connected to each other and therefore 
must be treated self-consistently to understand the complete picture 
of galactic chemical evolution. Semi-analytic models (SAMs) of 
galaxy formation follow cosmological halo growth and halo merg¬ 
ers and include physically and/or empirically motivated prescrip¬ 
tions of heating and cooling, star formation, metal enrichment, gas 


et al.|2006 
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Yates et al. 

20131 Lu et al.|201 1 |20141 Henriques et al.|20131 

2015 

Lu et al.||2015b|. They are much less computationally expensive 


to run than hydrodynamic simulations and are able to reproduce a 
number of galaxy properties for a broad range of stellar mass. How¬ 
ever, one major challenge for SAMs is simultaneously reproducing 
observed stellar masses, star formation rates (SFRs), and metallic- 
ities. The metallicities of low-mass galaxies are particularly sensi¬ 
tive to the galactic wind model because strong outflows are required 
to suppress star formation in low-mass systems (see e.g. |Lu et al.| 
|2014| for a detailed comparison and discussion). Moreover, even 
though different SAMs have been succcessfully tuned to match the 
z = 0 stellar mass function (SMF), many of them fail to match 
the observed the SMFs at high redshifts (e.g. |Somerville & Da^ 
|2014| >. At the same time, these models typically fail to match high- 
redshift MZR measurements and also diverge from one another in 
their MZR predictions. Nonetheless, it is encouraging that recently 
improved SAMs are able to reconcile stellar masses, colours, and 
SFRs of galaxies from z = 0-3 (e.g. |Henriques et al.|2013|[2015^ . 

Large-volume cosmological hydrodynamic simulations pro¬ 
duce large samples of galaxies and are powerful tools for statis¬ 


tical studies of galaxy properties (e.g. Bertone et al.|2007{|Dave et| 
|al.|2011b[|Torrey et al.||2014{|Schaye et al.|2015| . These simula¬ 

tions however usually have relatively poor mass and spatial resolu¬ 
tion. They cannot explicitly resolve the multi-phase structure of the 
interstellar medium (ISM), when and where star formation takes 
place, how galactic winds are launched by stellar feedback, and 
how the winds interact with the circum-galactic medium (CGM). 
Approximate, empirical “sub-grid” models of the ISM structure, 
star formation, and stellar feedback are required and used. For ex- 
ample, [Dave et al.| ( |201 la^ implemented a momentum-driven wind 
model, with wind mass loading factors and velocities prescribed as 
a function of bulk galaxy properties. In their implementation, hy¬ 
drodynamic interactions are temporarily suppressed as gas from the 
ISM is “kicked” into the galactic wind. Simulations using such sim¬ 
ple prescriptions reveal similar problems to the SAMs. [Torrey et af] 
( |2014^ found a steeper slope than observed at the low-mass end of 
the MZR. These authors attributed it to the low metal retention ef¬ 
ficiency in the presence of strong outflows, which were required in 
their model in order to prevent low-mass galaxies from forming too 
many stars. They further emphasized the tension between suppress¬ 
ing star formation and retaining enough metals in low-mass galax¬ 
ies. Furthermore, the star formation histories in these simulations 
are very different and not all consistent with observations at high 
redshifts. Many cosmological simulations tend to form too many 
stars at early times (e.g. |Dave et al.|2011a[ S parre et al.|2015[|Fiac-| 
|coni et al.|2015[ for a review, see | Somerville & Dave|2014t . Such 
problems are also common in SAMs. They are likely the result of 
imperfect star formation and stellar feedback models implemented 
in those simulations (cf. [Scannapieco et al.|2012) . Consequently, 
these simulations predict very different evolution of the MZR. 

Therefore, when using cosmological hydrodynamic simula¬ 
tions to understand the MZR and its evolution, one is required 
to capture the “correct” behavior of star formation, stellar feed¬ 
back, gas outflows, and the mixing and interaction of galactic winds 
with the CGM on all relevant scales. Encouragingly, [Obrej a et af] 
( |2014^ presented a suite of cosmological zoom-in simulations from 
the MaGICC project using an improved star formation and SNe 
feedback model. Their model includes an empirical prescription to 
approximate the effects of stellar feedback mechanisms operating 
before the first SNe explode. These authors showed that their simu¬ 
lations match the stellar mass-halo mass relation and the observed 
MZR from z = 0-3, for the eight galaxies in their sample. In this 
work, the first of a series, we will study the chemical evolution of 
galaxies using the FIRE (Feedb ack in Realistic Environment) sim¬ 
ulations I Hopkins et al. 2014) . The EIRE projecj^ is a series of 
cosmological zoom-in simulations that are able to follow galaxy 
merger history, interactions of galaxies with the IGM, and many 
other important processes. These simulations include a full set of 
realistic physical models and explicitly resolve the multi-phase 
structure of the ISM, star formation, and stellar feedback, with 
no need to tune parameters. The FIRE simulations successfully 
reproduce many observed galaxy properties, including the stellar 
mass-halo mass relation, star formation histories, the Kennicutt- 
Schmidt law, and the star-forming main sequence, from z = 0-6, for 


a broad range of galaxy masses in M, = 10‘*-10*’ Mq (Hopkins et 


|al.|2014| ). Also, the FIRE simulations predict reasonable covering 
fractions of neutral hydrogen in the halos of z = 2-3 Lyman Break 
Galaxies (LBGs; [Faucher-Giguere et al.|2015^ and self-consistently 
generate galactic winds with velocities and mass loading factors 
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broadly consistent with observational requirements ( [Muratov et al.| 
|2015| (. These results further justify the reliability to study galactic 
chemical evolution using the FIRE simulations. 

This paper focuses on the galaxy mass-metallicity relation. In 
companion papers, we will also study the stellar metallicity dis¬ 
tribution functions and [a/Fe] abundance ratio variation in dwarf 
galaxies, metallicity gradients and their origins, metal outflows and 
recycling. We start by describing the simulations in Section and 
present the mass-metallicity relation at different redshifts in Sec¬ 
tion]^ In Section]^ we discuss the key processes that drive the 
shape and evolution of the MZR. We summarize and conclude in 
SectionjS] 

2 THE SIMULATIONS 
2.1 Simulation Details 

This work is part of the FIRE project. A full description of the 
numerical methods and physics included in our simulations is pre¬ 
sented in ( [Hopkins et al.|2014| and references therein). We sum¬ 
marized their main features here. All the simulations use the newly 
developed GIZMO code ( [Hopkins[2014t in “P-SPH” mode. P-SPH 
adopts a Lagrangian pressure-entropy formulation of the smoothed 
particle hydrodynamics (SPH) equations ( [Hopkins][2013[ , which 
eliminates the major differences between SPH, moving-mesh, and 
grid codes, and resolves many well-known issues in traditional 
density-based SPH formulations. The gravity solver is a heavily 
modified version of the GADGET-3 code [Springel|2005[ l; and P- 
SPH also includes substantial improvements in the artificial vis¬ 
cosity, entropy diffusion, adaptive time-stepping, smoothing kernel, 
and gravitational softening algorithm. 

We use the multi-scale “zoom-in” initial conditions generated 
with the MUSIC code ( [Hahn & Abel[20I l) , using second-order La¬ 
grangian perturbation theory. The first set of simulations have been 
run down to z = 0 and cover halo masses IO^-IO'^Mq and stellar 
masses 10‘*-10*' Mq at z = 0 (luxx series). Most of them have been 
presented in [Hopkins et al.[ ( |2014[ . The simulations m09 and mIO 
are isolated dwarfs, constructed using the method from[ Onorbe et[ 
|^[2014| |. Simulations mIOv, mil, mllv, ml2q, ml2i, and ml3 
are chosen to match the initial conditions from the AGORA project 
( [Kim et ^[2014[ (, which will enable future comparisons with a 
wide range of simulation codes and physics implementations. Sim¬ 
ulation ml2v is based on the initial conditions studied in IKeres ^ 
[HemquistI ( [2009[ ( and [Eaucher-Giguere & Ker^ ( [2011[ . The sim¬ 
ulations with a label ‘v’ have relatively violent merger histories 
at z < 2, while the rest have more typical merger histories. The 
resolution of these simulations is chosen to scale with the mass 
of the system to ensure we are able to resolve the giant molecular 
clouds (GMCs). We also include a separate set of simulations run to 
z = 2 (z2hxxx series), which are presented in jFaucher-Giguere et[ 
[al.|2015> . Their main halos are chosen to host Lyman break galax¬ 
ies (LEG) and cover halo masses 1.9 x 10**-1.2 x atz = 2. 

Finally, we include another series of simulations only run to z ~ 6, 
but with extremely high resolutions (zSmxx series). These simula¬ 
tions are presented in [Ma et ^j2015^ . Their main halos cover halo 
masses from 7.7 x 10“-5.6 x lO'^M© at z = 6 and these galax¬ 
ies are believed to contribute most to the cosmic reionization (e.g. 
[Kuhlen & Faucher-Giguere 20 12[ [Robertson et al.pols} . The ini¬ 
tial conditions of all the simulations are summarized in Table[T] 

In our simulations, gas follows an ionized-(-atomic-(-molecular 
cooling curve from 10-10*° K, including metallicity-dependent 


fine-structure and molecular cooling at low temperatures and high- 
temperature metal-line cooling followed species-by-species for 11 
separately tracked species (H, He, C, N, O, Ne, Mg, Si, S, Ca, 
and Fe; see [Wiersma et alT][2009a^ . At each timestep, the ioniza¬ 
tion states and cooling rates are determined from a compilation of 
CLOUDY runs, including a uniform but redshift-dependent photo- 
ionizing background tabulated in [Faucher-Giguere et al.[ ( [2009[ l, 
and photo-ionizing and photo-electric heating from local sources. 
Gas self-shielding is accounted for with a local Jeans-length ap¬ 
proximation, which is consistent with the radiative transfer calcu¬ 
lations in [Faucher-Giguere et al.[p010[ l. 

Star formation is allowed only in dense, molecular, and self- 
gravitating regions with hydr ogen number density above some 
threshold nui = 10-100 cm“^ jHopkins et al.|[^13b| . Stars form 
at 100% efficiency per free-fall time when the gas meets these cri¬ 
teria. The self-gravity criterion is physically required to obtain the 
correct spatial star formation distribution in galaxies ([Hopkins et| 
[al.[2013b][Padoan & Nordlund|2011[ , but the galaxy-averaged star 
formation efficiency is regulated by feedback at much lower values 
(~ 1% per dynamical time, e.g. [Faucher-Giguere et al.|2013[ l. We 
stress that changing these parameters in a reasonable range only 
yields small and random variations to the global star formation his¬ 
tory, as long as feedback is active (see [Hopkins et al.[2011[[2012[ . 

Once a star forms, it inherits the metallicity of each tracked 
species from its parent gas particle. Every star particle is treated as 
a single stellar population with known mass, age, and metallicity. 
Then all the feedback quantities, including ionizing photon bud¬ 
gets, luminosities, stellar spectra, supemovae (SNe) rates, mechan¬ 
ical luminosities of stellar winds, metal yields, etc., are directly tab¬ 
ulated from the stellar population models in STAREURST99 ( [Lei-[ 
[therer et al.[1999^ , assuming a [Kroupa] ( [2002[ initial mass function 
(IMF) from O.I-lOO M(^ We account for several different stel¬ 
lar feedback mechanisms, including: (1) local and long-range mo¬ 
mentum flux from radiative pressure; (2) energy, momentum, mass 
and metal injection from SNe and stellar winds; and (3) photo¬ 
ionization and photo-electric heating. We follow [Wiersma et al.[ 
( [2009b[ l and include the metal yields from Type-II SNe, Type-I 
SNe, and stellar winds. We note that the Type-II SNe yield table 
from [Woosley & Weaver[ ( |1995[ adopted in our simulations pro¬ 
duce Mg roughly ~ 0.4 dex below the typical values in modem 
models (e.g. [Nomoto et al.[[2006l l. This will have little effect on 
the galaxy properties in our simulations, as Mg is not an important 
coolant. Nevertheless, we will add 0.4 dex to the Mg abundance to 
correct this in the analysis below. 

All simulations adopt a standard fiat ACDM cosmol¬ 
ogy with cosmological parameters consistent with Hq = 

70.2 km s“‘ Mpc"*, Da = 0.728, D„, = 1 - Da = 0.272, D* = 
0.0455, (Tg = 0.807 and n — 0.961 (e.g.[Hinshaw et al.|2013||Planck[ 
[Collaboration et al.[2013[ l. 

2.2 Halo Identification, Stellar Mass and Metallicity 

We use the Amiga Halo Finder (AHF; [Gill et al.|2004l[Knollmann[ 
[& Knebe[2009[ l to identify galactic halos and galaxies in our simu¬ 
lations. The AHF code uses the adaptive mesh refinement method 

^ In principle, the “IMF-averaged” approximation does not hold for the 
ultra-high resolution simulations in the zSmxx series, where the mass of 
a star particle is only 10-100 Mq. Nevertheless, we confirmed that these 
simulations predict similar global galaxy properties to those of much poorer 
resolutions (see[Ma et al.[2015^. 
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Table 1. Simulation Initial Conditions. 


Name 

^halo 

mt 



^dm 

Merger 

Notes 


(Me) 

(Me) 

(pc) 

(Me) 

(pc) 

History 


m09 

2.5e9 

2.6e2 

1.4 

1.3e3 

30 

normal 

isolated dwarf 

mlO 

O.SelO 

2.6e2 

3.0 

1.3e3 

30 

normal 

isolated dwarf 

mlOIr 

O.SelO 

2.1e3 

2.1 

1.0e4 

35 

normal 

low resolution 

mlOv 

O.SelO 

2.1e3 

7.0 

1.0e4 

70 

violent 

- 

mil 

1.4ell 

7.0e3 

7.0 

3.5e4 

70 

quiescent 

- 

mllv 

3.3ell 

5.6e4 

7.0 

3.0e5 

140 

violent 

- 

ml2v 

6.3ell 

3.9e4 

10 

2.0e5 

140 

violent 

several z < 2 mergers 

ml2q 

1.2el2 

7.1e3 

10 

2.8e5 

140 

late merger 

- 

ml2i 

l.lel2 

5.0e4 

14 

2.8e5 

140 

normal 

large (~ lORvir) box 

ml3 

6.0el2 

3.6e5 

21 

2.2e6 

210 

normal 

small group mass 

z2h3S0 

7.9ell 

5.9e4 

9 

2.9e5 

143 

normal 

- 

z2h400 

7.9ell 

5.9e4 

9 

2.9e5 

143 

quiescent 

- 

z2h450 

8.7ell 

5.9e4 

9 

2.9e5 

143 

normal 

- 

z2h506 

1.2el2 

5.9e4 

9 

2.9e5 

143 

violent 

- 

z2h550 

1.9ell 

5.9e4 

9 

2.9e5 

143 

quiescent 

- 

z2h600 

6.7ell 

5.9e4 

9 

2.9e5 

143 

violent 

- 

z2h650 

4.0ell 

5.9e4 

9 

2.9e5 

143 

normal 

- 

z2h830 

5.4ell 

5.9e4 

9 

2.9e5 

143 

normal 

- 

zSm09 

7.6e8 

16.8 

0.14 

81.9 

5.6 

quiescent 

ultra-high resolution 

zSmlO 

1.3el0 

131.6 

0.4 

655.6 

7 

normal 

ultra-high resolution 

zSmlOmr 

1.4el0 

l.le3 

1.9 

5.2e3 

14 

normal 

- 

zSmll 

5.6elO 

2.1e3 

4.2 

1.0e4 

14 

normal 

- 


Parameters describing the initial conditions for our simulations (units are physical): 

(1) Name: Simulation designation. 

(2) Mhaio- Approximate mass of the main halo at z = 0 (mxx series), z = 2 (z2hxxx series), or z = 6 
(zSmxx series). 

(3) Initial baryonic (gas and star) particle mass in the high-resolution region. 

(4) Minimum baryonic force softening (minimum SPH smoothing lengths are comparable or 
smaller. Force softening is adaptive (mass resolution is fixed). 

(5) Dark matter particle mass in the high-resolution region. 

(6) ^dm- Minimum dark matter force softening (fixed in physical units at all redshifts). 



MJM^ 

Figure 1. Stellar mass-gas-phase oxygen abundance relation at z = 0. The 
red solid and dashed curves represent the median and la dispersion of the 
SDSS MZR at z ~ 0.1 jTremonti et al.|2Q04^. T he open circles denote the 
data of the dwarf galaxy sample from |Lee et 5r]j2Q06) . Our simulations are 
in good agreement with observations fromM* = 10^-10^^ 3^0- 


and identifies halos and subhalos as groups of bound particles (dark 
matter, gas, and stars). In this work, we only consider those “well- 
resolved” halos that include more than 10^ bound particles, have at 
most 10% of their mass contaminated by low-resolution particles. 



Figure 2. Stellar mass-stellar metallicity relation at z = 0. The red solid 
and dashed curves are the median and \a dispersion of the SDSS MZR 
in the local universe JGallazzi et al.|2005|. The open circles represent the 
values of [Fe/H] of individual dwarfs from |Kirby et al.H2013) . Again, the 
agreement is good from IC^-IO^* 3^©- 

and contain at least 100 gas and 100 star particles, respectively. 
These criteria are somewhat arbitrary; but varying these numbers 
within a reasonable range will have little effect on our conclusions. 
If none of the halos meets these criteria in a snapshot (this happens 
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Figure 3. Stellar mass-gas-phase metallicity relation at all redshifts. Cyan dotted lines show the best linear lit log(Zgas/ZQ) = 12-|-log(0/H) — 9.0 = 
7 j[log(M,/MQ) — 10] -|-Zg 10 . The red dotted lines show the best fit for a fixed slope ■yg = 0.35. Note that a constant slope provides a very good fit, where 
the zero point evolves by ~ 1 dex from z = 0-6. 


in some snapshots at high redshifts (z ~ 6), where the galaxy pro¬ 
genitors are too small to contain so many particles), we will take 
the most massive halo in the high-resolution region in our analy¬ 
sis. We do not include subhalos/satellite galaxies in this work. The 
centre of a halo is located at the centre of mass of the finest refine¬ 
ment level. We adopt the virial overdensity from |Bryan & Norman| 
( |1998t , which evolves with cosmic time. 

We only consider the main galaxy in each halo. To remove the 
contamination of satellite galaxies, we exclude any gas/star parti¬ 
cle that is bound to a subhalo in the analysis below. We measure 
the galaxy stellar mass (M*) by summing over the mass of all star 
particles that belong to the main galaxy. Then we define its stel¬ 
lar metallicity (as well as the abundance of each tracked species) 
as mass-averaged metallicity (abundance) of all star particles. To 
separate halo gas and the ISM, we apply a simple temperature cri¬ 
teria and select all gas particles below 10"^ K as the ISM. In our 
simulations, this is equivalent to selecting gas above some density 
threshold of a few 0.1 cm“^ (we explicitly check the gas distribu¬ 
tion in the density-temperature plane), which is comparable to the 
mean gas density within a few stellar effective radii. It naturally 
picks warm ionized and cold neutral gas. We define the gas-phase 
metallicity as the mass-weighted metallicity of all gas particles that 


belong to the ISM (we compare and discu ss th ree different defini¬ 
tions of gas-phase metallicity in Appendix 


of all heavy elements in gas and stars, respectively. In Section 
we will primarily use oxygen abundance 12-|-log(0/H) to present 
gas-phase metallicities, which is defined in terms of number ratio 
of oxygen to hydrogen atoms, in order to directly compare with ob¬ 
servations. For stellar metallicity, we will primarily use Z* in the 
rest of this work. In the literature, gas-phase metallicity and stellar 
metallicity are also sometimes referred as Zgas and iron abundance 
[Fe/H] (in solar units), respectively. For these reasons, we provide 
the conversion between these quantities for our simulated galaxies. 
We will show the calibration in Appendix [B] but directly give the 
results here. For a solar metallicity of 0.02 and a solar iron abun- 


^ In many cosmological simulations with “sub-grid” models, gas-phase 
metallicity is usually defined as star-formation-rate-averaged metallicity. 
However, our simulations explicitly resolve multi-phase ISM structures and 
include realistic models of star' formation and feedback. Individual gas par¬ 
ticles are very sensitive to local feedback processes. For these reasons, we 
do not apply SF-averaged gas-phase metallicity to our simulations. 
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Figure 4. Stellar mass-stellar metallicity relation at all redshifts. Cyan dotted lines show the best linear ht at each redshift log(Zt/Z 0 ) = [Fe/H] + 0.2 = 
7 *[log(M*/M 0 ) — 10] +Z*,io. The red dotted lines show the best fit for a fixed slope 7 * = 0.40. Again, the slope is approximately constant, while the 
normalization evolves by ~ 1 dex. 


dance 0.00173 (both in mass fraction), we obtain 12 + log(0/H) = 
log(Zgas/Z0) +9.00 and [Fe/H] = log(Z*/Z0) -0.20. We empha¬ 
size that these relations may suffer from systematic uncertainties 
that originate from: (1) Type-II and Type-I SNe rates, (2) metal 
yields of tracked species from different channels, and (3) the solar 
abundances we adopt in our simulations. However, the shape and 
evolution of the MZR should be robust to these uncertainties. 


3 THE MASS-METALLICITY RELATION 

In this section, we present both the gas-phase and stellar MZR from 
z =0-6 and compare our results with observations and other sim¬ 
ulations. We will further explore the most important factors that 
shape the MZR and drive its evolution in the Sectionj^ 

3.1 The MZR at z = 0 

We begin by showing the gas-phase MZR at z = 0. In Figure 
we present the stellar mass-gas-phase oxygen abundance relation 
for our mxx series simulations at z = 0. For comparison, we also 
present the median and 2cr dispersion of the SDSS MZR from 


[Tremonti et al.| ( [2004| red solid and dashed lines) and the data of 
individual local dwarf galaxies compiled in |Lee et al.| ( |2006l open 
circles) in Figure[^ We remind the reader that these observed gas- 
phase oxygen abundances are derived from the relative strength of 
strong nebulae emission lines produced by photo-ionization from 
young massive stars, so that the observed gas-phase MZR only 
holds for star-forming galaxies. Also, we emphasize that the over¬ 
all shape of gas-phase MZR strongly depends on which empirical 
calibration it uses and the normalization of this relation differs by 
at most 0.7 dex between different calibrations ^Kewley & Ellison] 
|2008[ see also Figure]^. For these reasons, we do not apply any 
correction to these observed data but keep them in their original 
forms. 

In general, our simulations agree reasonably well with obser¬ 
vations across stellar mass from M* = 10®-10** Mq. However, our 
simulations do not show evidence for flattening at the high-mass 
end of the gas-phase MZR. The gas-phase metallicity increases 
with stellar mass up to M* ~ 10** Mq in our sample. The simula¬ 
tions predict^ligh^ higher metallicities than the observed relation 


from |Tremonti et al.| ( |2004| l atM* = 10 Mq. The most significant 
discrepancy is due to our ml3 run, which is a somewhat lower res¬ 
olution simulation of a massive galaxy and which did not include 
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Figure 5. The gas-phase and stellar metallicity at M* = 10*® Mq, Zj,io 
and Z* 10 a function of redshift. The solid lines are the best fit 
of exponential functions Zj io = 0.93exp(—0.43z) — 1.05 and Z, lo = 
0.67exp(-0.50z)-1.04. 


the possible effects of AGN feedback. For example, as it has been 
shown in |F[opkins et al.H2014^ , the main galaxy in ml3 have the 
cooling flow problem and never quenches at low redshift. The SFR 
of this galaxy is 3 Mq yr“*, which is fairly low in its star formation 
history, but significantly higher than observationally inferred values 
below z ~ 1. Consequently, this galaxy might be over-enriched at 
low redshift. If so, this suggests that additional physics, such as 
AGN feedback, is probably required to fully understand the chemi¬ 
cal evolution in massive galaxies, at least in the sense of quenching 
star formation. Alternatively, it has also been proposed that the ob¬ 
served MZR could continue to rise at the high-metallicity end when 
using new metallicity diagnostics that account for non-equilibrium 
electron energy distributions (see e.g. |Dopita et al.|2013[|Nicholls| 
|et al.|[20T^ . Furthermore, we note that the “flatness” of MZR at 
the high-mass end behaves very differently when applying differ¬ 
ent empirical calibrations (e.g. |Kewley & Ellison|2008) ). Therefore, 
we do not further quantitatively discuss the discrepancy between 
our simulations and observations at the massive-end of MZR, but 
rather focus on galaxies below M, = 10** Mq where our simula¬ 
tions are most robust. A larger sample of simulations with improved 
resolution at the massive end is required to make a robust compar¬ 
ison. 


Most of our simulated galaxies are still forming stars (at least 
very weakly) at z = 0, except for m09. The m09 is a low-mass 
isolated dwarf galaxy (comparable to the ultra faint dwarfs around 
the Milky Way), in which star formation has been shut down since 
z = 3 by cosmic reionization ( [Onorbe et al.|2015) l. At z = 0, this 
galaxy has lost almost all metals it produced (see Section]^. Al¬ 
though its gas-phase metallicity is an order of magnitude lower than 
the extrapolation of the observed MZR down to M* = 10“* Mq, 
it is not contradictory to observations in the sense that the gas- 
phase metallicity of such galaxies cannot be measured due to lack 
of strong nebular emission lines. 

In Figure]^ we show the stellar mass-stellar metallicity rela¬ 
tionship at z = 0 and compare our simulations with the SDSS sam¬ 
ple from |Gallazzi et aL |([2005[ red solid and dashed curves) and the 
dwarf galaxies from Kirby et al. 1 201^open circles). Note that the 
stellar metallicities from|Gallazzi et ah] j2005[> are measured from 


absorption features of galaxy-integrated spectra (mostly Mg and Fe 
lines), while the metallicities from |Kirby et al.H2013^ are derived 
from Fe abundances of individual stars. The conversion between 
different methods and their systematic uncertainties are complex 
and beyond the scope of this paper. For our purposes, we avoid any 
correction to these observations but present them in their original 
valuefl 

Our simulations match these observations quite well over the 
whole mass coverage from M* = 10‘*-10** Mq. The simulated 
sample shows a flatness in the stellar MZR around M* = 10** Mq 
at z = 0, consistent with the observed SDSS MZR from |Gallazzi| 
|et al.| ( |^05| l. This is the consequence of the fact that the growth 
of the more massive galaxies in our simulations is dominated by 
mergers and accretion of low-mass metal-poor satellites rather than 
in situ star formation at low redshifts. Therefore, the average stellar 
metallicities do not strongly evolve despite the fact that the stellar 
masses may grow considerably at low redshifts (see also|Choi et al.| 
|20T4l l. 


3.2 Evolution of the MZR 

Figure]^ and 1^ show the gas-phase and stellar MZR, respectively, 
from z = 0-6. We note that for z > 2 and z = 6 , we include the 
z2hxxx and zSmxx simulations in our analysis. The stellar MZR is 
tighter than the gas-phase MZR, i.e., the gas-phase MZR has larger 
scatter than stellar MZR at fixed stellar mass. This is because in 
our simulations, especially at high redshifts, star formation is domi¬ 
nated by multiple bursts, which drives bursts of gas outflows ( |Mura-| 
|tov et al.|2015[ ). As a consequence, instantaneous gas-phase metal¬ 
licities may have considerable time fluctuations associated with gas 
inflows, outflows, and mergers. This effect is larger at high red¬ 
shifts when the galaxy progenitors are of much lower masses and 
galaxy mergers are more common, resulting in some outliers that 
deviate from the main MZR at high redshifts. Despite the short- 
time-scale fluctuations, both the gas-phase and stellar metallicities 
increase with time on cosmological time-scales. At all times, gas- 
phase metallicities are higher than stellar metallicities, since gas- 
phase metallicities represent the current state of metal enrichment 
in the galaxies, while stellar metallicities reflect the average galac¬ 
tic metallicities across the whole time. Both metallicities should 
converge at high redshifts. 

To illustrate this quantitively, we fit the gas-phase and stel¬ 
lar MZR at different redshifts for our simulated galaxies with 
simple linear functions log(Zgas /Zq) = 12 + log(O/H)-9.0 = 
7 g[log(M*/M 0 ) - 10] +Zg,io and log(Z*/ZQ) = [Fe/H] +0.2 = 
7 ,[log(Mt/MQ) — 10] +Z*,io, where 7 j and 7 , are the slopes and 
Zg.io and Z,,io represent the typical gas-phase metallicity and stellar 
metallicity at M* = 10*® Mq. Although simple linear function do 
not capture the flatness of stellar metallicity above M* ~ 10** Mq 
at z < 1, it is sufficient for our purposes here. We use least-squares 
fitting to obtain the best fit (the cyan dotted lines in Figure]^ and 
1^. In principle, both the slopes and zero points should be func¬ 
tions of redshift. Nevertheless, the MZR at different redshifts have 
very similar slopes. For simplicity, we pick the mean slope of 
each relation and redo the linear fit using fixed slopes. We choose 
7 g = 0.35 and 7 * = 0.40 (red dotted lines in Figure and 


'* In Figure]^ we plot the values of [Fe/H] from |Kirby et al.|j2013) , avoid¬ 
ing the complicated conversion between [Fe/H] and Z*/Z 0 for the ob¬ 
served sample. 
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Figure 6. Stellar mass-gas-phase oxygen abundance relations at z = 0, 0.8, 2.2, and 3.0, as compared with a number of observations at these redshifts. In the 
upper panels, we show both the original relations (red lines) from|Tremonti et al.|^2004| z ~ 0) and|Zahid et al.|^2011| z ~ 0.8) and the relations converted to 
PP04 N2 calibration (cyan lines) following|Kewley & Ellison|^2008J. In the lower left panel, we show the observed MZR at z 2.3 from|Steidel et al.|^2014| 
the red line), [Sanders et al.|l|2015| the green line), and|Erb et al.|j2006| the yellow line). In the lower right panel, we show the best fitting from|Mannucci et| 
[^(2^z-3 .1). We also shift the |Erb et al.|j^Q6^ data downward by 0.4 dex for a comparison as motivated by Figure 5 in |Mannucci et ah|^2Q0^ . Our 
simulations are broadly consistent with observations over a wide range of stellar mass from z = 0-3, given the significant systematic uncertainties observational 
determinations of metallicities. 


and confirm that both the best linear fit and the fixed-slope fit de¬ 
scribe the simulations reasonably well. We then attribute the evo¬ 
lution of MZR to the evolution of Z^.io and Z,,io with redshift, 
which we show in Figure We fit these parameters as a func¬ 
tion of redshift by an exponential function F{z) = Ae,x\>[—Bz)+C. 
The best fit gives Z^jo = 0.93exp(—0.43z) — 1.05 and Z.^io = 
0.67 exp(—0.50z) — 1.04, respectively (the green and magenta lines 
in Figure]^. These give the gas-phase and stellar MZR from z = 0- 
6 as log(Zga.,/Z0) = 12 + log(0/H) - 9.0 = 0.35 [log(M./M0) - 
10] -h0.93exp(-0.43z) - 1.05 and log(Z./Z0) = [Fe/H] + 0.2 = 
0.40 [log(M*/M0) — 101 + 0.67 exp(—0.50z) — 1.04, respectively. 


In general, the fitting functions above represent the gas-phase 
and stellar MZR fairly well for our simulated galaxies, except for 
the flattening of the stellar MZR above M* ~ 10*' Mq at z = 0. 
We emphasize that these results have systematic uncertainties from 
Type-II and Type-la SNe rates, the solar abundance, and the metal 
yield tables we implement in our simulations. When using these 
fitting functions, one should notice the uncertainties and make ad¬ 
justments accordingly. 


3.3 Comparison with Observations and Other Models 

In Figure]^ we compare the gas-phase MZR between our simula¬ 
tions and a number of observations at multiple redshifts. We show 
the observed MZR at z ~ 0 jTremonti et al.|2004t , z ~ 0.8 Zahid 
et al.|20lT| l, z ~ 2.2 jSteidel et al.|2014 [Sanders et al.|2015| Erb et 


al.|2006t , and z ~ 3.1 ( jMannucci et al. 


2009 | l. We recall that these 


observed relations are originally obtained using different calibra¬ 
tions and the systematic uncertainty between different metallicity 
diagnostics could be up to 0.7 dex ^Kewley & Ellison||2008| ). To 
illustrate this point, we also convert all the observed relation to the 
N2 calibration from jPettini & Pag^p004| PP04 hereafter) unless 
their original data are already presented using this calibration. For 
[Tremonti et al.| ( [2004t and [Zahid et al.| ( [201 1[ >, we do the conver¬ 
sion following the formula from Kewley & Ellison| ( [2008[ Table 3 
therein). In either case, we present both their original relations and 
the converted relations using PP04 N2 calibration in Figure At 
z ~ 2.2, the observed relations are at already presented in PP04 N2 
calibration (e.g. Steidel et al. 2014[[Sanders et al.|2015[[Erb et al.[ 


[2006[ l. [Mannucci et al.[ ( [2009 i adopted a very different metallicity 
calibration, which is established using z > 3 galaxy samples only. 
Figure 5 in [Mannucci et al.]([2009[ suggests that the MZR evolves 
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Figure 7. Stellar mass-gas-phase oxygen abundance relation at z = 0, 0.8, 2.2, and 3.0, as compared with other numerical simulations and semi-analytic 
models. We renormalize other works to 12 -|- log(0/H) = 8.9 at M* = IO'^Mq at z = 0 with respect to our simulations. Red and green lines show the results 
from cosmological simulations presented in |Torrey et al.| j20T4j and |Dave et al.| J201 lb|, respectively, which used popular “sub-grid” models for galactic 
winds. Magenta, cyan, and yellow lines show the predictions of three semi-analytic models from |Lu et al.|j20T4] the Lu model, the Somerville model, and 
the Croton model, respectively). All of these models reproduce the correct z = 0 stellar mass function, but none of them correctly reproduces the slope or the 
redshift evolution of the MZR. 


by ~ 0.4 dex from z ~ 3.1 to z ~ 2.2. Motivated by their results, 
we also move the z ~ 2.2 MZR from |Erb et al.|j?006t downward 
by 0.4 dex for a comparison (lower right panel in Figure]^. 

In general, our simulations are in reasonable agreement with 
these observations in a broad range of stellar mass at z = 0-3, es¬ 
pecially when the observed relations are in their original forms. We 
emphasize that the empirical calibrations developed from the local 
universe are not necessarily valid for high-redshift galaxies (e.g. 
ISteidel et al.|2014[|Kewley et al.|2015) ). Given the large systematic 
uncertainties, we do not provide a detailed quantitative discussion 
of the discrepancies between our simulations and observations. Our 
results on the evolution of the MZR in Section |3^ are predictions 
that can be tested more accurately as our understanding of the ob¬ 
servational systematic uncertainties improves. 

In Figure [7] we also compare the MZR from our simulations 
with other cosmological simulations and semi-analytic models. We 
compare our results with two other simulations, |Torre y et al.|j2014| 
red lines) and |Dave et ^r]j201 1b[ g reen lines), and three semi- 
analytic models from Lu et al. ( 2014[ the Lu model, magenta; the 
Somerville model, cyan; the Croton model, yellow). These mod¬ 
els adopt “sub-grid” empirical models of galactic winds and stellar 
feedback, which couple some fraction of energy and/or momen¬ 
tum from SNe to the gas, and force certain amount of the gas to 


escape the galaxy. Note that the metal yields and solar abundance 
used in different works are not exactly the same, we renormalize 
all the z = 0 MZR to 12 -f log(0/H) = 8.9 at M, = 10'° Mq for 
comparison. At z = 0, |Torrey et al.| ) (2014| l and the Lu model show 
Steeper slopes at the low-mass end, due to the low metal retention 
efficiency in low-mass galaxies, a consequence of invoking strong 
outflows to suppress star formation in these galaxie^ Some models 
predict higher metallicities at the most massive end. Furthermore, 
these models show significant discrepancies at z > 2. Our simula¬ 
tions predict much stronger evolution of MZR from z = 3-0 than 
any other models. Particularly, the Somerville model and the Cro- 


^ This can be simply illustrated using the “leaky box” model (e.g.|Schmidt| 
|1963) . Assuming the outflow rate is proportional to the star formation rate 
(Afoul = f] • SFR, where 77 is the mass loading factor), the metallicity is in¬ 
versely proportional to 1 -h 77 . Low-mass galaxies are very efficient in driv¬ 
ing outflows and thus have high mass loading factors compared to massive 
galaxies. In SAMs and some simulations with “sub-grid” feedback models, 
it is often assumed that either the metals are well mixed in the system or 
that the outflowing gas has a metallicity comparable to the metallicity in the 
ISM. As a consequence, low-mass galaxies tend to lose a large fraction not 
only of their gas but also of their metals, and therefore to end up with very 
low metallicities. 
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Figure 8 . Metal retention fraction for our simulated galaxies at z = 0. Mz (< 
f?vir) is the total amount of metals retained (in both gas and stars) within the 
virial radius. yM* (y is the mean effective yield) is the total metal mass 
produced by stars in the galaxies. The retained fraction of metal in the halo 
increases with stellar mass, from 30% at M* = 10® Mq to about unity at 
M, > 10'° Mq. However, the ultra-faint dwarfs (e.g. m09) are only able to 
retain 2% of their metals in the halo. 


ton model predict inverse evolution trends - the gas-phase metal- 
licity decreases at lower redshifts at fixed stellar mass - in contrast 
with observations and other models. We recall that although these 
models are tuned to match the observed stellar mass function at 
z = 0, they tend to predict systematically higher stellar mass func¬ 
tions than the observed ones for M* < 10*' Mq at z > 0 
|& Dave|2014 ^, a consequence of the fact that galaxies in these mod¬ 
els form too many stars at early time (e.g. |Dave et al.|201 la[ |Sp a^ 
|et al.||2015| |Fiacconi et al.||2015) . In Section we further ex¬ 
plore how the different star formation histories between these mod¬ 
els cause the discrepancies in the MZR at high redshifts. 


Somerville 


4 DISCUSSION 

We showed above that the gas-phase and stellar MZR in our sim¬ 
ulations agree broadly with available observations at different red¬ 
shifts. We also found that our predictions diverge significantly from 
those of several large-volume cosmological hydrodynamical simu¬ 
lations and semi-analytic models. In this section, we explore the 
key factors that drive the shape and evolution of the MZR and dis¬ 
cuss why our predictions differ from some other models. 


4.1 Where are the metals? 


Our simulations produce much higher metallicities for galaxies of 


stellar mass M* < 10 Mq than Torrey et al. 


12014b and the Lu 


model in |Lu et ah] ( |2014| l, indicating that our low-mass galaxies re¬ 
tain more metals compared to those models, despite the fact that 
these galaxies have high wind mass loading factors up to 100. To 
explicitly show this, we present in Figure the metal mass frac¬ 
tion retained within Rvir as a function of stellar mass for the simu¬ 
lated sample at z = 0. The numbers are obtained as follows. First, 
we estimate the effective yield y for every simulation as the ratio 
between total metal mass (in both gas and stars) and the total stel¬ 
lar mass in the whole simulation volume. Then the metal retention 


fraction for a galaxy is simply the ratio between the total metal 
mass within the virial radius, Mz{< Rvir), and yMt, where M* is 
the galaxy stellar mass as defined in Section [2^ Thus, yM, rep¬ 
resents the total amount of metal ever produced by the stars in the 
galaxy. As shown in Figure]^ the metal retention fraction generally 
increases with stellar mass. In our simulated sample, galaxies above 
Mt = 10*** ® Mq are able to keep almost all metals they have pro¬ 
duced. At much lower masses (M* = 10®-10’ Mq), they can still 
retain at least 30% to a half of their metals within the halo. In con¬ 
trast, the ultra-faint dwarf in our sample, m09 (M* = 4 x 10** Mq), 
only retains 2% of its metals within Rvir at z = 0. 

To quantify in more detail how metals are retained in galaxy 
halos, we show in Figure the cumulative metal retention frac¬ 
tion, as a function of radius, for different gas phases (cool gas with 
T <10* K and warm gas with 10'* K < T < 4 x 10® K:^ At z = 0 
(top row), low-mass galaxies such as mlO (M, = 2 x 10° Mq) have 
most of their metals in the warm CGM, while in massive galax¬ 
ies like mI2i (M* = 6 x 10*** Mq), the majority of the metals are 
found in stars. This trend is qualitatively consistent with the em¬ 
pirical halo metal budget presented in |Peeples et al.|p014| Fig. 6). 
In most cases, we find that only a small fraction of the total metal 
mass is found in hotter (T > 4 x 10® K) gas. Our results are in con¬ 
trast with the large-volume simulations of |Ford et al.| ( |2015| l based 
on a parameterized galactic wind model, in which stars, ISM, and 
the cool CGM contain comparable metal masses for halos of mass 
similar to our ml2i run. 

Regarding the spatial distribution of metals, in mil (M* = 
2x10® Mq), over 60% of the metals are concentrated in the central 
0.1 Rvir (mostly the galaxy) and only a small fraction (< 40%) of 
metals are in the circumgalactic medium (CGM) or lost into the 
IGM. In massive systems such as ml2i (M, = 6 x 10*** Mq), almost 
all the metals are in the central 0.1 Rvir. In low-mass galaxies like 
mlO (M* = 2 X 10® Mq), metals are more evenly distributed among 
the galaxy, the CGM, and the IGM. In ultra-faint dwarfs like m09 
(Mf = 4 X 10'* Mq ), most of the metals it has ever produced are lost 
in the IGM by z = 0. This is consistent with the fact that outflows in 
low-mass galaxies are more efficient (they have much higher mass 
loading factor) and can propagate more easily to large radii than 
in massive systems jMuratov et al.||^15b . It also shows that the 
metals are far from well mixed in the halo of more massive galaxies 
with stellar mass M, > 10® Mq (or in terms of halo mass M„ir > 
10“ Mq). 

For a comparison, we also show the cumulative metal distribu¬ 
tion for the progenitors of these galaxies at z = 3 (the bottom panel 
in Figure]^. Similar to z = 0, a significant fraction of metals are 
still retained in Rvir at z = 3, although metals are more uniformly 
distributed from the centre to a few virial radii. These galaxies have 
much lower mass than their low-redshift decedents, and thus they 
are more efficient in driving gas outflows from star-forming regions 
throughout the halo. 


4.2 Circum-galactic O VI 

Although this paper is primarily focused on the metallicity of gas 
and stars inside galaxies, it is useful to check whether our sim¬ 
ulations are consistent with observed CGM metal absorption. In 


® In our simulations, most of the diffuse (n^{ <0.1 cm“^) gas has tempera¬ 
ture T > 10'^ K so a temperature cut at T = lO'^ K also effectively separates 
ISM and CGM gas, justifying our approach of using gas with T < 10“* K to 
evaluate gas-phase ISM metallicities. 
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Figure 9. Cumulative metal mass in selected simulated halos at z = 0 (top) and z = 3 (bottom), normalized by the total metal mass produced by stars (yM*). 
The red dashed, blue dotted, orange dash-dotted, and black solid lines show the metal mass in stars, cool gas (T < 10** K), warm gas (10** K < T < 4 X 10^ 
K), and total, respectively. The stellar mass of each galaxy is indicated at the top left corner of each panel and the black dotted lines show the virial radius. At 
z = 0, most of the metals in our more massive simulated galaxies such as ml2i are in stars and within 0.1 i?vir of the halo center, while in low-mass galaxies, 
the majority of metals are found in the warm CGM. In low-mass galaxies at z = 0 and in high-redshift galaxies, a larger fraction of the metals are found at 
larger radii from the halo center, consistent with the fact that galactic outflows are more powerful in these systems. 



kpc 

Figure 10. 0 VI column density map for the ml2i halo at z = 0. We crudely 
assume that a fraction /ovi = 0.2 of the oxygen is in O VI and only include 
warm and hot gas (T > 10“* K) in the halo. The characteristic Vovi drops 
from ~ 10*^ cm~^ at impact parameter b = 20 kpc from the central galaxy 
to ~ 10*^'^ cm~^ at b = 200 kpc. The simulation agrees well with the 
O VI columns measured by COS-Halos [Tumlinson et aklpOll) around 
low-redshift ~ L* star-forming galaxies at impact parameters b < 50 kpc 
but appeal's to underestimate O VI columns by a factor of a few at larger 
impact parameters. Overall the agreement with observed O VI columns is 
reasonable given the uncertainties in ionization correction. 


addition to the overall metal budget discussed above, the COS- 
Halos program has provided useful measurements of O VI ab¬ 
sorption around ~ L* galaxies at z « 0.1 — 0.4 jTumlinson et al.| 
|2011| >. Figure [T0| shows the O VI column density map around our 
ml2i simulated halo at z = 0. For this comparison, we assume that 
a fraction /ovi = 0.2 of the oxygen is in O VI and only include 
warm and hot gas {T > 10"* K) in the halo, /ovi = 0.2 is the maxi¬ 
mum expected if the oxygen is in collisional ionization equilibrium, 
though it is possible that OVI is also photoionized and/or subject 
to non-equilibrium effects (e.g., |Oppenheimer & Dave|200^ |Op-| 
Ipenheimer & Schaye||2013) so that this ionization fraction is not 
a strict upper limit. The figure shows that for ml2i the character¬ 
istic Now drops from ~ 10*^ cm“^ at impact parameter b = 20 
kpc from the central galaxy to ~ 10*^'^ cm“^ at b = 200 kpc. The 
simulation agrees well with the O VI columns measured by |Tum-| 
llinson et al.| ( |2011) around low-redshift ~ L* star-forming galaxies 
at impact parameters i < 50 kpc but appears to underestimate O VI 
columns by a factor of a few at larger impact parameters. Overall 
the agreement with observed O VI columns is reasonable given the 
uncertainties in ionization correction. More systematic and detailed 
comparisons of CGM metal statistics from the FIRE simulations 
with observations will be reported in future papers (Hafen et al., in 
preparation). 


4.3 Metal outflows, inflows, and recycfing 

SAMs and large-volume cosmological simulations require “sub¬ 
grid” models of galactic winds, which often incorporate fairly crude 
approximations. In this subsection, we further examine the metal 
inflow and outflow rates and the metallicities of gas inflows and 
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Time (Gyr) Time (Gyr) 

Figure 11. Upper. Metal inflow (blue) and outflow rates (red) from z = 0—4. Solid and dotted lines show the metal inflow/outflow rates measured at 0.25 i?vir 
and Rvi[., respectively. Bottom'. Metallicities of inflowing/outflowing gas. The black line shows the metallicity of the ISM. All quantities are averaged over a 
time-scale of 400 Myr. Metals are efficiently ejected in fountains reaching 0.25 Rvir. but they do not usually reach Ryir - they are either deposited in the halo 
or recycled efficiently in galactic fountains. Outflowing gas that escapes from the halo at Ry;, tends to be less enriched than the gas in the ISM. 
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Figure 12. Left: Stellar mass fraction /* = M*/(Mgas +M*) as a function of stellar mass. Mgas here is the total gas mass in the halo (not only in the galaxy). 
Middle'. Stellar metallicity Z* as a function of /*. Right: Gas-phase metallicity Zgas as a function of /*. For consistency, Zgas here is the average metallicity of 
all gas in the halo (including both the ISM and the halo gas). Black points and red points show the primary FIRE simulations at z = 0 and z = 3, respectively. 
Blue dotted lines show the simple “closed box” model predictions assuming an effective metal yield of y = 0.02. The z = 0 and z = 3 galaxies share the same 
Z*-/* and Zgas-/* relations, but the relation evolves by ~ 0.5 dex from z = 3-0. This indicates that the evolution of the MZR is associated with the 

evolution of /* (at a fixed stellar mass) at different redshifts. The major offset between our simulations and the predictions of the “closed box” model is largely 
due to the fact that the metals are not perfectly mixed throughout the halo. Especially in massive galaxies, gas tends to be more metal-enriched in the central 
star-forming regions than in the outer halo, so stellar metallicities tend to be higher and gas-phase metallicities (including the halo gas) are lower than the 
predictions of the “closed box” model. 
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outflows in our simulations and compare with the assumptions of 
common “sub-grid” models. 

We follow [Faucher-Giguere et j20 1 1"^ and [Muratov et alT] 
( [2015^ and define the gas outflow rates, metal outflow rates, and 
metallicities of outflow gas as 


dM ^ r , 

i ' ' 

(1) 

dMmefiX -> ^ r/ lAi 

Q, -E-|,|Z.M/dL, 

/ ' ' 

(2) 

^outflow r\ / ? 

at at 

(3) 


where M, and Z; are the mass and metallicity of the i* gas parti¬ 
cle within the shell of thickness dL = 0.1 f^vii with radial velocity 
outwards if- > 0. The inflow rates and inflow metallicities are 
defined in the same way but for gas particles with inward radial 
velocity v*- ^ < 0. The upper panels in Figure 11 show the metal 
inflow/outflow rates at 0.25 (blue/red solionnes) and at 
(blue/red dotted lines) for our mlO (left) and mil (right) simula¬ 
tions. We average the inflow/outflow rates on a time-scale of 400 
Myr. In either case, the net metal outflow rates are considerably 
lower at l^vir than at 0.25 IJvir, indicating that the metals are either 
deposited in the halo or returned back to the ISM. At high red- 
shifts, metals ejected in outflows can be more easily driven to Ry„ 
than at low redshifts. At 0.25 R^ir, metal inflow rates are compa¬ 
rable to metal outflow rates, suggesting a high efficiency of metal 
recycling. The lower panels in Figure [TT] show the average metal¬ 
licities of inflows and outflows at both 0.25 IJvir and at l^vir, as com¬ 
pared to the metallicity of the ISM (black solid lines). The outflow 
metallicities are much lower at f?vir than at 0.25 Rm (and even more 
so than in the ISM), because outflowing gas sweeps up and mixes 
with more metal-poor gas in the halo as it propagates outwards. 
This is particularly important for low-mass galaxies, such as mlO 
(M, = 2 X 10^ Mq), which can have wind mass loading factors up 
to ~ 100, yet retain a large fraction of the metals they produced in 
their halos. 

Our analysis calls into question a number of assumptions and 
approximations often adopted in analytic, semi-analytic, and large- 
volume cosmological hydrodynamic models of galaxy formation. 
First of all, unlike often assumed in analytic and semi-analytic 
models, metals are generally not well-mixed in galaxy halos (e.g. 
Figure]^. In particular, in many “sub-grid” galactic wind models, 
wind gas is assumed to have a metallicity directly related to the ISM 
metallicity (e.g. [Dave et al.|2011a[|Torrey et al.|2014^ , an assump¬ 
tion that oversimplifies the complex mass and metal loading that 
takes places in our more explicit simulations. Our simulations also 
indicate that metal re-accretion onto galaxies (recycling) is impor¬ 
tant on small scales, an effect which is not well captured in semi- 
analytic models and in “sub-grid” models that either assume that 
the ejected gas never returns to the galaxy, or which ignore hydro- 
dynamical interactions between the wind and the gas close to the 
galaxy. 

Recently, |Lu et al!] ( |2015b^ compared three different SAM 
feedback models — one including only gas ejection, one including 
both gas ejection and recycling, and the other including a model of 
“preventive” feedback. |Lu et al.| j2015b| l found that none of these 
models could simultaneously reproduce the MZR, the distribution 
of metals in different phases, and the SFR observed at z = 0-3. 
This finding is consistent with the picture suggested by our high- 
resolution simulations that the chemical evolution of galaxies is a 


complex process and that it is necessary to self-consistently model 
galaxy-halo interactions in order to capture it faithfully. It is en¬ 
couraging that our cosmological simulations with explicit stellar 
feedback and hydrodynamical interactions tracked at all times ap¬ 
pear to produce a low-mass-end slope of the MZR that is closer 
to observations than most previous models, without the need for 
parameter tuning. Our results are broadly consistent with those of 
[Brook et al.| ( |2014} , who also highlighted the importance of metal 
mixing with the COM and recycling for explaining the MZR. The 
simulations of [Brook et al. j20M l also provide a fair match to the 


observed MZR at z = 0-3 ( Obreja et al.|2014^ ). 


4.4 Why the MZR evolves with redshift? 


Another major difference between our simulations and other the¬ 
oretical work is we predict much stronger evolution of the MZR 
from z = 3-0 (e.g. the stellar metallicity increases by 0.5 dex at 
fixed stellar mass, see Figure]^. Observations and some theoretical 
models suggest a fundamental metallicity relation (FMR) between 
stellar mass, star formation rate, and metallicity that holds for star¬ 
forming galaxies both in the local universe and at high redshifts 


(e.g. 

Mannucci et al. |2010| 2Ull||Lilly et al.||2013||Obreja et al.| 

|20T4 

Cullen et al.|[2014 Zahid et al.||2014|l. Motivated by these 


results, we attempt to qualitatively illustrate what might be the pri¬ 
mary factor that drives the evolution of MZR in this section. We 
start by reviewing the simplest chemical evolution model, i.e., the 
“closed box” model, which predicts the stellar and gas-phase metal¬ 
licities as a function of stellar mass fraction, /* = Mt/(Mgas -l-M*) 
as the following 


Z, 


1 -/* 

/* 


ln(l-/*) + ! 


-y in(i-/.), 


(4) 

(5) 


where y is the effective metal yield (e.g. [ Schmidt|1963[ [Talbo t &] 
[Amett|1971[|Searle & Sargent|1972l ). The parameter /* describes 
the fraction of baryons that have been turned into stars, and 1 — /* 
is the “gas fraction”. In Figure [T^ we show the relation between 
stellar and gas-phase metallicities and /*, respectively (the mid¬ 
dle and right panels), for our mxx series simulations at z = 0 and 
z = 3 (black and red points). We emphasize that we account for both 
the halo gas and the ISM in the total gas mass when calculating 
ft, since halo gas is actively involved in supplying star formation 
and metal exchange in most cases. For consistency, the gas-phase 
metallicities shown in Figure[T^is the average metallicity of all gas 
in the halo. For illustrative purpose, we also show the simple pre¬ 
dictions from the “closed box” model, assuming an effective metal 
yield of y = 0.02 (blue dotted lines in Figure [T^. 

The simulated data at z = 0 and z = 3 overlap with each other 
in the Z*-/, and Zgas-/* diagrams. In the left panel of Figure \T2\ 
we also show the relation between /* and M, for these galaxies at 
both redshifts. There is a systematic offset (~0.5 dex) in the /*-M* 
relation between galaxies at z = 0 and 3. Note that in the limit of 
ft 1, one has Z,,Zgas oc /,. Therefore, the 0.5 dex offset in /*- 
Mt relation propagates to the 0.5 dex evolution of the MZR from 
z = 3 to 0. This suggests that the evolution of the MZR is associated 
with the evolution of /* (at a fixed stellar mass) within the halo at 
different redshifts, providing a first hint of a universal metallicity 
relation between stellar mass, gas mass, and metallicities (cf. [Both-'| 
[well et ar|[2013[ [Zahid et ^[2014[ for observational evidences). 
In simulations with “sub-grid” feedback models and semi-analytic 
models, where the z = 0 stellar mass functions are tuned to match 
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observations, galaxies tend to form a large fraction of their stars at 
high redshift and therefore their evolution is weaker at lower red- 
shift (e.g. [Somerville & Da^|2014^ , as opposed to observations 
and our simulations. In other words, these models produce higher 
/* than our simulations at fixed stellar mass at z > 0 and an /*-M* 
relation barely evolving from z = 3-0. Therefore, galaxies in those 
models are more metal-enriched at high redshifts and the evolution 
of the MZR is weaker than our simulations. 

Our simulations are qualitatively consistent with the simple 
“closed box” predictions applied to halo quantitie^ This is not 
unreasonable because a large fraction (order unity) of metals are 
retained within the virial radius at both redshifts (see e.g.. Figure 
1^. However, we emphasize that one should not think our simulated 
galaxies are closed boxes, because the metals are not perfectly well- 
mixed in the galactic halo. This explains the major offset between 
the “closed box” model and our simulations (middle and right pan¬ 
els in Figure [T^, especially in the most massive systems where 
this effect is stronger. Since gas in the centre of the galaxy tends to 
be more metal-enriched than gas in the outer halo and stars pref¬ 
erentially form in the central region, stellar metallicities tend to be 
higher and the gas-phase metallicities (including the halo gas) are 
lower than the predictions of the closed box model (applied to halo 
quantities). The mixing of metals is very complex and associated 
with galactic fountains on different scales. Although the “closed 
box” model gives a natural relation between stellar mass, gas mass, 
and the metallicities, the parameterization of a universal metallic- 
ity relation for galactic quantities (i.e., excluding the halo) is more 
complicated than the simple model. This is worth further investiga¬ 
tion in more detail in future work. 


5 CONCLUSION 

We use a series of high-resolution cosmological zoom-in simula¬ 
tions spanning halo masses 10^-10*^ Mq and stellar masses lO"*- 
10*' Mq at z = 0 from the FIRE project to study the galaxy mass- 
metallicity relations at z = 0-6. These simulations include explicit 
models of multi-phase interstellar medium, star formation, and stel¬ 
lar feedback. As has been shown in previous papers, these sim¬ 
ulations successfully reproduce many observed galaxy properties, 
including the stellar mass-halo mass relation, star-forming main 
sequence, the Kennicutt-Schmidt law, star formation histories, etc., 
for a wide range of galaxies at many redshifts ( [Hopkins et al.|2014l >. 
These simulations also predict reasonable covering fractions of 
neutral hydrogen in the halos of z = 2-3 LBGs ( [Faucher-Giguere ^ 
|al.[2015^ and self-consistently generate galactic winds with veloci¬ 
ties and mass loading factors broadly consistent with observational 
requirements ( [Muratov et al.[20r5) . These simulations adopt “stan¬ 
dard” stellar population models and metal yield tables from Type-I 
and Type-II supernovae and stellar winds, following species-by- 
species for 11 separately tracked elements. Our key conclusions 
include the following. 

(i) The simulations predict galaxy mass-metallicity relations 
that agree reasonably well with a number of observations from z = 
0-3 for a broad range of stellar masses. Both gas-phase and stellar 

^ We emphasize that in the analog of Figure[^where we measure /* using 
only the gas in the galaxy (i.e., excluding the halo gas), all the galaxies are 
well below the predictions of the closed box model and there is no well- 
defined relation, indicating that galaxies themselves are far from closed 
boxes. This suggests the necessity of accounting for halo gas as reservoirs 
in galaxy evolution. 


metallicities evolve monotonically from z = 0-6, with higher metal 
abundance at low redshifts at fixed stellar mass. The best linear fits 
of the MZR for our simulated galaxies as a function of redshift 
are log(Zgas/Z0) = 12 + log(0/H) - 9.0 = 0.35 [log(M*/M©)- 
10]-F0.93exp(-0.43z)- 1.05 and log(Z./Z0) = [Fe/H] + 0.2 = 
0.40 [log(Mt/Af0) — 10]+0.67exp(—0.50z)— 1.04, for gas-phase 
metallicity and stellar metallicity, respectively. We emphasize that 
the normalizations may have systematic uncertainties that originate 
from the SNe rates, yield tables, and solar abundance we adopt, but 
the evolution of the MZR is robust to these uncertainties. 

(ii) The stellar MZR becomes flat around M, ~ 10*' Mq since 
z = 0, because the most massive galaxies in our simulations evolve 
via mergers and accretion of satellites rather than in situ star for¬ 
mation at low redshifts. Therefore, the stellar metallicity does not 
increase despite the fact that the stellar mass grows considerably. 
We do not see the flatness in the gas-phase MZR at the high-mass 
end seen in observations, because gas continues to be enriched by 
non-negligible star formation. This apparent discrepancy may be 
due to the more limited resolution in our ml3 run or to the lack 
of AGN feedback in our simulations. AGN might be required to 
quench star formation below z ~ 1 in such massive galaxies. 

(iii) The evolution of MZR is associated with the evolution of 
the gas/stellar mass fraction within the inner halo (not just inside 
the galaxy effective radius) at different redshifts. This provides a 
first hint of a universal metallicity relation between stellar mass, 
gas mass, and metallicities, but its parameterization for galactic 
quantities (as opposed to for halo quantities, which behave more 
like a closed box) is much more complicated than simple analytic 
models. We will investigate this in more detail in future work. 

(iv) Galaxies above 10** Mq can retain a large fraction of their 
metals in the halo even up to z = 3. The net metal outflow rates 
near the virial radius are always lower than those near the galaxy, 
indicating that the metals either get deposited in the halo or return 
back to the ISM. The high metal inflow rates and the high metallic¬ 
ity of inflowing gas at 0.25 Rvir suggest a high efficiency of metal 
recycling (a finding that we have confirmed using particle tracking; 
Angles-Alcazar et ah, in prep.). On average, the outflows at outer 
radii are much less metal-enriched than those at the inner radius. 
This effect helps resolve the tension between the need for strong 
gas outflows and high metal retention fractions in low-mass galax¬ 
ies. 

(v) These differential recycling and metal retention effects are 
not properly accounted for in most semi-analytic and early gen¬ 
eration of “sub-grid” feedback models that are popular in cosmo¬ 
logical simulations. As a result, these simplified models cannot si¬ 
multaneously reproduce the galaxy mass function and the slope and 
redshift evolution of the MZR. By explicitly resolving the “missing 
physics” in these models, we reconcile the long-standing discrep¬ 
ancy, and provide a clear way forward to improve the sub-grid and 
semi-analytic models. 

Nevertheless, our simulations are still limited in sample size. 
In the near future, we will expand our simulations to include more 
dwarf galaxies covering halo mass fromMhaio = 10*-10'* Mq and 
to enlarge our sample at the most massive end to better understand 
whether the flattening of the MZR is real and what drives the flat¬ 
ness. This may depend critically on AGN feedback. We will pro¬ 
vide quantitative analysis on metal outflow rates, outflow metallic¬ 
ities, metal recycling, and their relation with galaxy properties in 
future work (Muratov et ah, in preparation; Angles-Alcazar et al, 
in preparation). 
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APPENDIX A: DIFFERENT DEFINITIONS OF 
GAS-PHASE METALLICITY 

In this work, the gas-phase metallicity is defined as the mass- 
weighted average metallicity of all gas particles below lO"* K, 
which we refer as the ISM gas. In principle, there are many alterna¬ 
tive approaches to define gas-phase metallicities. In this section, we 
discuss three definitions and compare them with each other: (1) the 
average metallicity of all gas particles below lO’* K in the galaxy 
(our default definition), (2) the average metallicity of all gas parti¬ 
cles within 0.1 Rvk, and (3) the average metallicity of all gas parti¬ 
cles with temperature between 7,000-15,000 K and density above 
0.5 cm“^. In Figure[^ we compare definition (1) and (2) in the left 
panel and (1) and (3) in the right panel for all galaxies presented in 
Figurej^ 

Definition (1) is designed to automatically select all the warm 
ionized gas and cold neutral gas (the ISM), definition (2) aims to 
pick the gas in the star-forming regions, and definition (3) is ob- 
servationally motivated to select the nebular gas which produce the 
strong nebular emission lines in star-forming galaxies. In general, 
these definitions are consistent with each other. Most of the galax¬ 
ies lie very close to the y — x relation in each panel of Figure 
However, there are a few outliers in these diagrams. Definition (2) 
can be problematic in merging systems, where the halo centre may 
deviate far from the stellar bulk and thus 0.1 RvL does not neces¬ 
sarily probe the star-forming region. Definition (3) is largely af¬ 
fected by abundance variance between gas particles, since there are 
usually not many gas particles at any single instant that meet the 
temperature and density criteria. However, a time-averaged version 
of definition (3) removes most of the outliers. Therefore, we argue 
that our default definition is more adaptive and flexible than other 
definitions. 


APPENDIX B: METALLICITIES IN DIFFERENT FORMS 

In this work, we primarily use 12-f log(0/H) and Z* to present 
gas-phase metallicity and stellar metallicity, respectively. In the lit¬ 
erature, gas-phase metallicity and stellar metallicity are sometimes 
presented in terms of ZgiiS and [Fe/H]. Therefore, we also provide 
the conversion between these different forms of metallicities for 
comparison. We emphasize these conversions are obtained from 
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Figure 1. Gas-phase oxygen abundances in different definitions. Left: The relation of gas oxygen abundances between definition (1) the average metallicity 
of all gas particles below 10"^ K and (2) the average metallicity of all gas particles within 0.1 Right: The relation of gas oxygen abundances between 
definition (1) and (3) the average metallicity of all gas particles with temperature between 7,000-15,000 K and density above 0.5 cm~^. The cyan dashed lines 
show the y = X relation. The black points show all the data presented in Figure]^ Different definitions agree well, and have no qualitative effect on any of our 
conclusions. Most of the “outliers” are caused by transient, stochastic time variability. 
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Figure 2. Relations between different forms of metallicities. Left: Gas-phase oxygen abundance 12 + log(0/H) vs. gas-phase metallicity (mass fraction of 
all metals) Zgas- Right: Stellar iron abundance [Fe/H] vs. stellar metallicity Z*. Black dots collect all the data points presented in this work. The cyan lines 
represent the best fits of these relations with slope unity. These definitions give essentially identical results, and are equivalent, for all of our results in this 
paper. 


our simulations only and there are systematic uncertainties origi¬ 
nating from the uncertain relative metal yields between species and 
solar abundances we adopt. 

In Figure]^ we show the relations between 12 + log(0/H) 
and log(Zgas/Z0) (left panel) and the relation between [Fe/H] and 
log(Z*/ZQ) (right panel), where we adopt a solar metallicity Z0 = 
0.02 and a solar iron abundance of 0.00173, both in mass fraction. 
In both panels, we collect data of all the simulated galaxies at all 
epochs we present earlier in this paper. Both relations are extremely 
tight and have slope unity, which ensures the validity, at least to 
the first order, to use either quantity to represent metallicities inter¬ 
changeably. The best fits for our simulations are 12-|-log(0/H) = 
log(Zga.s/Z0) +9.0 and [Fe/H] = log(Z*/Z0) - 0.20. We empha¬ 
size that there relations may suffer from systematic uncertainties 
that originate from: (1) Type-II and Type-I SNe rates, (2) metal 


yields of tracked species from different channels, and (3) the solar 
abundances we adopt in our simulations. 
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